# Slovakian elections

check.packages <- function(pkg){
  new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
  if (length(new.pkg))
    install.packages(new.pkg, dependencies = TRUE,
                     Repos = "https://cran.rstudio.com")
  sapply(pkg, require, character.only = TRUE)
}

pkgList <- c("tidyverse", "lme4","broom","ggeffects","rstudioapi",
             "texreg","arm","broom.mixed","ggthemes",
             "interplot", "corrplot")
check.packages(pkgList)


# 
setwd("~/ArcCR500/Slovensko")

library(readxl)
VOLBY_SVK <- read_excel("~/ArcCR500/Slovensko/Slovenské volby 2020.xlsx")
View(VOLBY_SVK)


# Korelační matice

KOR_MATICE <- VOLBY_SVK[, c("Ucast_20", "OLANO_20","PS_a_SPOLU_20",
                             "Kotlebovci_Ľudová_strana_Naše_Slovensko_20", "Sloboda_a_Solidarita_20", "SME_RODINA_20",
                             "SMER_20", "ZA_LUDI_20", "VZ_Vysokoskolske", "Vek_61", "EA_nezamestnaný", "RS_rozvedení",
                            "OLANO_16",	"SME_RODINA_16", "MOST_HÍD_16", "SNS_16", "SMER-SD_16", "KDH_16", "ĽS_Naše_Slovensko_16", "SIEŤ_16",	"SaS_16")]

COR_ZAKLADNI <-cor(KOR_MATICE, use="complete.obs", method="pearson")
colnames(COR_ZAKLADNI) <- c("Účast", "Olano", "PS/Spolu",  
                            "Kotleba", "SaS", "Sme Rodina", "SMER", "Za ĽUDÍ", "Vysokoškoláci",
                            "Důchodci", "Nezaměstnanost", "Rozvedení",
                            "Olano 16", "Sme Rodina 16", "Most-Híd 16", "SNS 16", "SMER 16",
                            "KDH 16","Kotleba 16", "Siet 16", "SaS 16")

rownames(COR_ZAKLADNI) <-  c("Účast", "Olano", "PS/Spolu",  
                             "Kotleba", "SaS", "Sme Rodina", "SMER", "Za ĽUDÍ", "Vysokoškoláci",
                             "Důchodci", "Nezaměstnanost", "Rozvedení",
                             "Olano 16", "Sme Rodina 16", "Most-Híd 16", "SNS 16", "SMER 16",
                             "KDH 16","Kotleba 16", "Siet 16", "SaS 16")
library(corrplot)
corrplot(COR_ZAKLADNI, method="number", type="upper", tl.col = "black",  number.cex = 0.5, cl.cex = 0.5, tl.cex = 0.5)

# Sociodemografie

KOR_MATICE <- VOLBY_SVK[, c("Ucast_20", "VZ_Vysokoskolske", "Vek_61", "EA_nezamestnaný", "RS_rozvedení", "VV_VERICI","NAR_maďarská","Vel_obyv_LN")]
COR_ZAKLADNI <-cor(KOR_MATICE, use="complete.obs", method="pearson")
corrplot(COR_ZAKLADNI, method="number", type="upper", tl.col = "black") #,  number.cex = 0.8, cl.cex = 0.8, tl.cex = 0.8)


# Regression models


PS  <- lm(PS_a_SPOLU_20 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + NAR_maďarská, data=VOLBY_SVK, weights=Zapsani_volici)
summary (PS)

library(coefplot)
coefplot(PS, intercept = FALSE,  title = "zp: PS/Spolu", xlab = "Nestandardizované hodnoty koeficientu  + 95% k.i.",
         ylab = "Název  koeficientu", sort = c("magnitude"), pointSize = 3, color = "purple", 
         newNames=c(VZ_Vysokoskolske="Vysokoškoláci (%)",
                    Vek_61="Důchodci (%)",
                    EA_nezamestnaný="Nezaměstnaní (%)",
                    RS_rozvedení="Rozvedení (%)",
                    NAR_maďarská="Podíl Maďarů (%)")) + 
  geom_vline(xintercept=0, linetype="dashed", color = "red", size=1.1)



OLANO  <- lm(OLANO_20 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + NAR_maďarská, data=VOLBY_SVK)
library(coefplot)
coefplot(OLANO, intercept = FALSE,  title = "zp: OLANO", xlab = "Nestandardizované hodnoty koeficientu  + 95% k.i.",
         ylab = "Název  koeficientu", sort = c("magnitude"), pointSize = 3, color = "yellow3", 
         newNames=c(VZ_Vysokoskolske="Vysokoškoláci (%)",
                    Vek_61="Důchodci (%)",
                    EA_nezamestnaný="Nezaměstnaní (%)",
                    RS_rozvedení="Rozvedení (%)",
                    NAR_maďarská="Podíl Maďarů (%)")) + 
  geom_vline(xintercept=0, linetype="dashed", color = "blue", size=1.1)



Kotleba  <- lm(Kotlebovci_Ľudová_strana_Naše_Slovensko_20 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + NAR_maďarská, data=VOLBY_SVK)
library(coefplot)
coefplot(Kotleba, intercept = FALSE,  title = "Zp: Kotleba", xlab = "Nestandardizované hodnoty koeficientu  + 95% k.i.",
         ylab = "Název  koeficientu", sort = c("magnitude"), pointSize = 3, color = "brown", 
         newNames=c(VZ_Vysokoskolske="Vysokoškoláci (%)",
                    Vek_61="Důchodci (%)",
                    EA_nezamestnaný="Nezaměstnaní (%)",
                    RS_rozvedení="Rozvedení (%)",
                    NAR_maďarská="Podíl Maďarů (%)")) + 
  geom_vline(xintercept=0, linetype="dashed", color = "red", size=1.1)


SME_RODINA  <- lm(SME_RODINA_20 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + NAR_maďarská, data=VOLBY_SVK)
library(coefplot)
coefplot(SME_RODINA, intercept = FALSE,  title = "Zp: SME RODINA", xlab = "Nestandardizované hodnoty koeficientu  + 95% k.i.",
         ylab = "Název  koeficientu", sort = c("magnitude"), pointSize = 3, color = "red", 
         newNames=c(VZ_Vysokoskolske="Vysokoškoláci (%)",
                    Vek_61="Důchodci (%)",
                    EA_nezamestnaný="Nezaměstnaní (%)",
                    RS_rozvedení="Rozvedení (%)",
                    NAR_maďarská="Podíl Maďarů (%)")) + 
  geom_vline(xintercept=0, linetype="dashed", color = "blue", size=1.1)




####################################################################################################################################
####################################################################################################################################
####################################################################################################################################
VOLBY_SVK <- read_excel("~/ArcCR500/Slovensko/Slovenské volby 2020.xlsx")


# OLS
####################################################################################################################################
PS  <- lm(PS_a_SPOLU_20 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská + Vel_obyv_LN, data=VOLBY_SVK, weights=Zapsani_volici, weight=Zapsani_volici)

OLANO  <- lm(OLANO_20 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská  + Vel_obyv_LN, data=VOLBY_SVK, weight=Zapsani_volici)

Kotleba  <- lm(Kotlebovci_Ľudová_strana_Naše_Slovensko_20 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská  + Vel_obyv_LN, data=VOLBY_SVK, weight=Zapsani_volici)

SME_RODINA  <- lm(SME_RODINA_20 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská  + Vel_obyv_LN, data=VOLBY_SVK, weight=Zapsani_volici)

PS_a_SPOLU <- lm(PS_a_SPOLU_20 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská  + Vel_obyv_LN, data=VOLBY_SVK, weight=Zapsani_volici)

Sloboda_a_Solidarita <- lm(Sloboda_a_Solidarita_20 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská  + Vel_obyv_LN, data=VOLBY_SVK, weight=Zapsani_volici)

Slovenská_národná_strana <- lm(Slovenská_národná_strana_20 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská  + Vel_obyv_LN, data=VOLBY_SVK, weight=Zapsani_volici)

SMER <- lm(SMER_20 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská  + Vel_obyv_LN, data=VOLBY_SVK, weight=Zapsani_volici)

ZA_LUDI <- lm(ZA_LUDI_20 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská  + Vel_obyv_LN, data=VOLBY_SVK, weight=Zapsani_volici)

MOSTHÍD <- lm(MOSTHID_20 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská  + Vel_obyv_LN, data=VOLBY_SVK, weight=Zapsani_volici)

KDH <- lm(Kresťanskodemokratické_hnutie_20 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská  + Vel_obyv_LN, data=VOLBY_SVK, weight=Zapsani_volici)



library(texreg)
htmlreg(list(OLANO, SMER, SME_RODINA, Kotleba, PS_a_SPOLU, Sloboda_a_Solidarita, ZA_LUDI, KDH, MOSTHÍD, Slovenská_národná_strana), 
        file = "./SVK_2020.html", stars = numeric(0),
        custom.model.names = c("Olano","SMER", "SME_RODINA","LSNS","PS SPOLU",
        "SaS",  "ZA LUDI", "KDH", "MOST_HID",
        "SNS"), 
        reorder.coef = c( 2, 3, 4, 5, 6, 7, 8, 1),
        custom.coef.names = c("(Constant)","College degree (%)" ,  "61+ years (%)","Unemployed (%)" , "Divorced", "Catholics (%)",
                              "Hungarians (%)", "Number of inhabitants (log)"), 
        custom.gof.names = c("R2","Adj.R2","Number of cases", "RMSE"), 
        digits = 3, # decimal places for coefficients and SEs
        single.row = FALSE) # display BETA and SEs under each other





# HLM
####################################################################################################################################
VOLBY_SVK <- read_excel("~/ArcCR500/Slovensko/Slovenské volby 2020.xlsx")
VOLBY_SVK$OKRES <- as.factor(VOLBY_SVK$OKRES)

library(lme4)

OLANO  <- lmer(OLANO_20 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská  + Vel_obyv_LN + (1| KOD_OKRESU), data=VOLBY_SVK  , weight=Zapsani_volici)
summary(OLANO)

PS  <- lmer(PS_a_SPOLU_20 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská + Vel_obyv_LN + (1 |KOD_OKRESU ), data=VOLBY_SVK  , weight=Zapsani_volici)

SMER <- lmer(SMER_20 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská  + Vel_obyv_LN  + (1| KOD_OKRESU), data=VOLBY_SVK, weight=Zapsani_volici)

SME_RODINA  <- lmer(SME_RODINA_20 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská  + Vel_obyv_LN  + (1| KOD_OKRESU), data=VOLBY_SVK, weight=Zapsani_volici)

Kotleba  <- lmer(Kotlebovci_Ľudová_strana_Naše_Slovensko_20 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská  + Vel_obyv_LN  + (1| KOD_OKRESU), data=VOLBY_SVK, weight=Zapsani_volici)

PS_a_SPOLU <- lmer(PS_a_SPOLU_20 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská  + Vel_obyv_LN  + (1| KOD_OKRESU), data=VOLBY_SVK, weight=Zapsani_volici)

Sloboda_a_Solidarita <- lmer(Sloboda_a_Solidarita_20 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská  + Vel_obyv_LN  + (1| KOD_OKRESU), data=VOLBY_SVK, weight=Zapsani_volici)


ZA_LUDI <- lmer(ZA_LUDI_20 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská  + Vel_obyv_LN+ (1| KOD_OKRESU), data=VOLBY_SVK  , weight=Zapsani_volici)

KDH <- lmer(Kresťanskodemokratické_hnutie_20 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská  + Vel_obyv_LN  + (1| KOD_OKRESU), data=VOLBY_SVK, weight=Zapsani_volici)

MOSTHÍD <- lmer(MOSTHID_20 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská  + Vel_obyv_LN + (1| KOD_OKRESU), data=VOLBY_SVK, weight=Zapsani_volici)

Slovenská_národná_strana <- lmer(Slovenská_národná_strana_20 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská  + Vel_obyv_LN  + (1| KOD_OKRESU), data=VOLBY_SVK, weight=Zapsani_volici)

library(texreg)
htmlreg(list(OLANO, SMER, SME_RODINA, Kotleba, PS_a_SPOLU, Sloboda_a_Solidarita, ZA_LUDI, KDH, MOSTHÍD, Slovenská_národná_strana), 
        file = "./SVK_2020_HLM.html", stars = numeric(0),
        custom.model.names = c("Olano","SMER", "SME_RODINA","LSNS","PS SPOLU",
                               "SaS",  "ZA LUDI", "KDH", "MOST_HID",
                               "SNS"), 
        reorder.coef = c( 2, 3, 4, 5, 6, 7, 8, 1),
        custom.coef.names = c("(Constant)","College degree (%)" ,  "61+ years (%)","Unemployed (%)" , "Divorced", "Catholics (%)",
                              "Hungarians (%)", "Number of inhabitants (log)"), 
        custom.gof.names = c("AIC","BIC","Log Likelihood", 
                     "Num. obs.", "Num. groups: District",
                     "Var: District (Intercept)","Var: Residual"), 
        digits = 3, # decimal places for coefficients and SEs
        single.row = FALSE) # display BETA and SEs under each other

OLANO  <- lmer(OLANO_20 ~ 1 + (1| KOD_OKRESU), data=VOLBY_SVK, weight=Zapsani_volici)

# Bez velikosti obce
OLANO  <- lmer(OLANO_20 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská  + (1| KOD_OKRESU), data=VOLBY_SVK  , weight=Zapsani_volici)
summary(OLANO)
18.36 / (18.36+39484.01)*100


PS  <- lmer(PS_a_SPOLU_20 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská + (1 |KOD_OKRESU ), data=VOLBY_SVK  , weight=Zapsani_volici)

SMER <- lmer(SMER_20 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská  + (1| KOD_OKRESU), data=VOLBY_SVK, weight=Zapsani_volici)

SME_RODINA  <- lmer(SME_RODINA_20 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská  + (1| KOD_OKRESU), data=VOLBY_SVK, weight=Zapsani_volici)

Kotleba  <- lmer(Kotlebovci_Ľudová_strana_Naše_Slovensko_20 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská  + (1| KOD_OKRESU), data=VOLBY_SVK, weight=Zapsani_volici)

PS_a_SPOLU <- lmer(PS_a_SPOLU_20 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská  + (1| KOD_OKRESU), data=VOLBY_SVK, weight=Zapsani_volici)

Sloboda_a_Solidarita <- lmer(Sloboda_a_Solidarita_20 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská  + (1| KOD_OKRESU), data=VOLBY_SVK, weight=Zapsani_volici)


ZA_LUDI <- lmer(ZA_LUDI_20 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská  + (1| KOD_OKRESU), data=VOLBY_SVK  , weight=Zapsani_volici)

KDH <- lmer(Kresťanskodemokratické_hnutie_20 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská  + (1| KOD_OKRESU), data=VOLBY_SVK, weight=Zapsani_volici)

MOSTHÍD <- lmer(MOSTHID_20 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská  +  (1| KOD_OKRESU), data=VOLBY_SVK, weight=Zapsani_volici)

Slovenská_národná_strana <- lmer(Slovenská_národná_strana_20 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská  + (1| KOD_OKRESU), data=VOLBY_SVK, weight=Zapsani_volici)

library(texreg)
htmlreg(list(OLANO, SMER, SME_RODINA, Kotleba, PS_a_SPOLU, Sloboda_a_Solidarita, ZA_LUDI, KDH, MOSTHÍD, Slovenská_národná_strana), 
        file = "./SVK_2020_HLM_BEZ_LNOBYV.html", stars = numeric(0),
        custom.model.names = c("Olano","SMER", "SME_RODINA","LSNS","PS SPOLU",
                               "SaS",  "ZA LUDI", "KDH", "MOST_HID",
                               "SNS"), 
        reorder.coef = c( 2, 3, 4, 5, 6, 7, 1),
        custom.coef.names = c("(Constant)","College degree (%)" ,  "61+ years (%)","Unemployed (%)" , "Divorced", "Catholics (%)",
                              "Hungarians (%)"),
        custom.gof.names = c("AIC","BIC","Log Likelihood", 
                             "Num. obs.", "Num. groups: District",
                             "Var: District (Intercept)","Var: Residual"), 
        digits = 3, # decimal places for coefficients and SEs
        single.row = FALSE) # display BETA and SEs under each other




# HLM 2016
####################################################################################################################################


OLANO  <- lmer(OLANO_16 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská  + Vel_obyv_LN + (1| KOD_OKRESU), data=VOLBY_SVK  , weight=Zapsani_volici)
summary(OLANO)

SMER <- lmer(SMER_SD_16 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská  + Vel_obyv_LN  + (1| KOD_OKRESU), data=VOLBY_SVK, weight=Zapsani_volici)

SME_RODINA  <- lmer(SME_RODINA_16 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská  + Vel_obyv_LN  + (1| KOD_OKRESU), data=VOLBY_SVK, weight=Zapsani_volici)

Kotleba  <- lmer(ĽS_Naše_Slovensko_16 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská  + Vel_obyv_LN  + (1| KOD_OKRESU), data=VOLBY_SVK, weight=Zapsani_volici)

Sloboda_a_Solidarita <- lmer(SaS_16 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská  + Vel_obyv_LN  + (1| KOD_OKRESU), data=VOLBY_SVK, weight=Zapsani_volici)

KDH <- lmer(KDH_16 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská  + Vel_obyv_LN  + (1| KOD_OKRESU), data=VOLBY_SVK, weight=Zapsani_volici)

MOSTHÍD <- lmer(MOST_HÍD_16 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská  + Vel_obyv_LN + (1| KOD_OKRESU), data=VOLBY_SVK, weight=Zapsani_volici)

Slovenská_národná_strana <- lmer(SNS_16 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská  + Vel_obyv_LN  + (1| KOD_OKRESU), data=VOLBY_SVK, weight=Zapsani_volici)

library(texreg)
htmlreg(list(OLANO, SMER, SME_RODINA, Kotleba, Sloboda_a_Solidarita, KDH, MOSTHÍD, Slovenská_národná_strana), 
        file = "./SVK_2016_HLM.html", stars = numeric(0),
        custom.model.names = c("Olano","SMER", "SME_RODINA","LSNS",
                               "SaS",  "KDH", "MOST_HID",
                               "SNS"), 
        reorder.coef = c( 2, 3, 4, 5, 6, 7, 8, 1),
        custom.coef.names = c("(Constant)","College degree (%)" ,  "61+ years (%)","Unemployed (%)" , "Divorced", "Catholics (%)",
                              "Hungarians (%)", "Number of inhabitants (log)"), 
        custom.gof.names = c("AIC","BIC","Log Likelihood", 
                             "Num. obs.", "Num. groups: District",
                             "Var: District (Intercept)","Var: Residual"), 
        digits = 3, # decimal places for coefficients and SEs
        single.row = FALSE) # display BETA and SEs under each other



#  Změna podílu hlasů
####################################################################################################################################




OLANO  <- lmer(OLANO_ROZDIL ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská  + Vel_obyv_LN + (1| KOD_OKRESU), data=VOLBY_SVK  , weight=Zapsani_volici)

SMER <- lmer(SMER_ROZDIL ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská  + Vel_obyv_LN  + (1| KOD_OKRESU), data=VOLBY_SVK, weight=Zapsani_volici)

Kotleba  <- lmer(KOTLEBA_ROZDIL ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská  + Vel_obyv_LN  + (1| KOD_OKRESU), data=VOLBY_SVK, weight=Zapsani_volici)

Sloboda_a_Solidarita <- lmer(SAS_ROZDIL ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská  + Vel_obyv_LN  + (1| KOD_OKRESU), data=VOLBY_SVK, weight=Zapsani_volici)

KDH <- lmer(KDH_ROZDIL ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská  + Vel_obyv_LN  + (1| KOD_OKRESU), data=VOLBY_SVK, weight=Zapsani_volici)

MOSTHÍD <- lmer(MOST_ROZDIL ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská  + Vel_obyv_LN + (1| KOD_OKRESU), data=VOLBY_SVK, weight=Zapsani_volici)

Slovenská_národná_strana <- lmer(SNS_ROZDIL ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská  + Vel_obyv_LN  + (1| KOD_OKRESU), data=VOLBY_SVK, weight=Zapsani_volici)

SME_RODINA  <- lmer(SME_RODINA_16 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + RS_rozvedení + VV_VERICI + NAR_maďarská  + Vel_obyv_LN  + (1| KOD_OKRESU), data=VOLBY_SVK, weight=Zapsani_volici)


library(texreg)
htmlreg(list(OLANO, SMER, SME_RODINA, Kotleba, Sloboda_a_Solidarita, KDH, MOSTHÍD, Slovenská_národná_strana), 
        file = "./SVK_2016_2020_DIFF_HLM.html", stars = numeric(0),
        custom.model.names = c("Olano","SMER", "SME_RODINA","LSNS",
                               "SaS",  "KDH", "MOST_HID",
                               "SNS"), 
        reorder.coef = c( 2, 3, 4, 5, 6, 7, 8, 1),
        custom.coef.names = c("(Constant)","College degree (%)" ,  "61+ years (%)","Unemployed (%)" , "Divorced", "Catholics (%)",
                              "Hungarians (%)", "Number of inhabitants (log)"), 
        custom.gof.names = c("AIC","BIC","Log Likelihood", 
                             "Num. obs.", "Num. groups: District",
                             "Var: District (Intercept)","Var: Residual"), 
        digits = 3, # decimal places for coefficients and SEs
        single.row = FALSE) # display BETA and SEs under each other





# Náhodný koeficient
####################################################################################################################################
####################################################################################################################################


library(broom)

#Rozvedení
OLANO  <- lmer(OLANO_20 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + 
                 RS_rozvedení + VV_VERICI + NAR_maďarská  + (1 + RS_rozvedení| NAZEV_OKRES), data=VOLBY_SVK  , weight=Zapsani_volici)
summary(OLANO)

# Koeficient
ranef(OLANO) %>%
  filter(variable=="RS_rozvedení") %>%
  dplyr::select(c(level, estimate, std.error)) %>%
  # Keep only columns identifying the country, the Beta, and SE
  rename(NAZEV_OKRES = level) %>%
  mutate(estimate = estimate + fixef(OLANO)["RS_rozvedení"]) %>%
  # Add fixed-effect to deviation from overall slope so as to get
  # the slope in each country
  ggplot(aes(x = reorder(NAZEV_OKRES, -estimate),
             y = estimate)) +
  geom_point(size = 2, color = "black") +
  labs(x = "District",
       y = "Effect of the divorced (%) on OLANO vote share") +
  geom_errorbar(aes(ymin = estimate - 1.64*std.error,
                    ymax = estimate + 1.64*std.error),
                size = 1,
                width = 0, color = "black") +
  geom_hline(yintercept = fixef(OLANO)["RS_rozvedení"],
             size = 2,
             linetype = "dashed",
             color = "black") +
  geom_hline(yintercept = fixef(OLANO)["RS_rozvedení"] +
               1.64*se.fixef(OLANO)["RS_rozvedení"],
             size = 1,
             linetype = "dotted") +
  geom_hline(yintercept = fixef(OLANO)["RS_rozvedení"] -
               1.64*se.fixef(OLANO)["RS_rozvedení"],
             size = 1,
             linetype = "dotted") + theme_base() +
  theme(axis.text.x = element_text(angle = 45, vjust = 0.95, hjust = 1)) +
  theme(axis.text.x = element_text(size = 10.5))




#Nezaměstnanost
OLANO  <- lmer(OLANO_20 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + 
                 RS_rozvedení + VV_VERICI + NAR_maďarská  + (1 + EA_nezamestnaný| NAZEV_OKRES), data=VOLBY_SVK  , weight=Zapsani_volici)
summary(OLANO)

# Koeficient
ranef(OLANO) %>%
  augment(ci.level = 0.90) %>%
  filter(variable=="EA_nezamestnaný") %>%
  dplyr::select(c(level, estimate, std.error)) %>%
  # Keep only columns identifying the country, the Beta, and SE
  rename(NAZEV_OKRES = level) %>%
  mutate(estimate = estimate + fixef(OLANO)["EA_nezamestnaný"]) %>%
  # Add fixed-effect to deviation from overall slope so as to get
  # the slope in each country
  ggplot(aes(x = reorder(NAZEV_OKRES, -estimate),
             y = estimate)) +
  geom_point(size = 2, color = "black") +
  labs(x = "District",
       y = "Effect of unemployment (%) on OLANO vote share") +
  geom_errorbar(aes(ymin = estimate - 1.64*std.error,
                    ymax = estimate + 1.64*std.error),
                size = 1,
                width = 0, color = "black") +
  geom_hline(yintercept = fixef(OLANO)["EA_nezamestnaný"],
             size = 2,
             linetype = "dashed",
             color = "black") +
  geom_hline(yintercept = fixef(OLANO)["EA_nezamestnaný"] +
               1.64*se.fixef(OLANO)["EA_nezamestnaný"],
             size = 1,
             linetype = "dotted") +
  geom_hline(yintercept = fixef(OLANO)["EA_nezamestnaný"] -
               1.64*se.fixef(OLANO)["EA_nezamestnaný"],
             size = 1,
             linetype = "dotted") + theme_base() +
  theme(axis.text.x = element_text(angle = 45, vjust = 0.95, hjust = 1)) +
  theme(axis.text.x = element_text(size = 10.5))






# "OLANO_20","PS_a_SPOLU_20", "Kotlebovci_Ľudová_strana_Naše_Slovensko_20", "Sloboda_a_Solidarita_20", "SME_RODINA_20", "SMER_20"


#Nezaměstnanost
OLANO  <- lmer(Kotlebovci_Ľudová_strana_Naše_Slovensko_20 ~ VZ_Vysokoskolske + Vek_61 + EA_nezamestnaný + 
                 RS_rozvedení + VV_VERICI + NAR_maďarská  + (1 + EA_nezamestnaný| NAZEV_OKRES), data=VOLBY_SVK  , weight=Zapsani_volici)
summary(OLANO)

# Koeficient
ranef(OLANO) %>%
  broom::augment(ci.level = 0.90) %>%
  dplyr::filter(variable=="EA_nezamestnaný") %>%
  dplyr::select(c(level, estimate, std.error)) %>%
  # Keep only columns identifying the country, the Beta, and SE
  rename(NAZEV_OKRES = level) %>%
  mutate(estimate = estimate + fixef(OLANO)["EA_nezamestnaný"]) %>%
  # Add fixed-effect to deviation from overall slope so as to get
  # the slope in each country
  ggplot(aes(x = reorder(NAZEV_OKRES, -estimate),
             y = estimate)) +
  geom_point(size = 2, color = "black") +
  labs(x = "District",
       y = "Effect of unemployment (%) on SaS vote share") +
  geom_errorbar(aes(ymin = estimate - 1.64*std.error,
                    ymax = estimate + 1.64*std.error),
                size = 1,
                width = 0, color = "black") +
  geom_hline(yintercept = fixef(OLANO)["EA_nezamestnaný"],
             size = 2,
             linetype = "dashed",
             color = "black") +
  geom_hline(yintercept = fixef(OLANO)["EA_nezamestnaný"] +
               1.64*se.fixef(OLANO)["EA_nezamestnaný"],
             size = 1,
             linetype = "dotted") +
  geom_hline(yintercept = fixef(OLANO)["EA_nezamestnaný"] -
               1.64*se.fixef(OLANO)["EA_nezamestnaný"],
             size = 1,
             linetype = "dotted") + theme_base() +
  theme(axis.text.x = element_text(angle = 45, vjust = 0.95, hjust = 1)) +
  theme(axis.text.x = element_text(size = 10.5))
